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A least squares technique has been applied to periodic earth temperature data for the purpose 
of determining basic characteristics of earth temperature cycles, such as thermal diffusivity, average 
temperature, amplitude, and phase angle of the temperature cycle. A new procedure was developed 
for obtaining a single thermal diffusivity which represents an average over time and depth at a particu- 
lar temperature site. This thermal diffusivity was obtained as a nonlinear part of least squares con- 
stants which yielded a best-fit harmonic curve to a given set of observed earth temperatures. The 
thermal diffusivity thus calculated and the calculation method developed are preferred to those obtained 
by current practice, which yields two thermal diffusivities, one based on amplitude decay and another 
on phase angle shift. 
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Nomenclature 

Aq, Ai, At, . . . A N = Least squares constants or best 
unbiased estimates of «o, «i, «2, . . . «v (°C) 

a () , ai, «2, . . • as — Cosine coefficients of expected 
earth temperature function (°C) 

Z?i, #2, B,u • • • fi.v = Least squares constants or best 
unbiased estimates of fei, 62, 63, ■ • • bx (°C) 

61, b 2 , 63, • • • />.\ = Sine coefficient of expected 

A earth temperature function (°C) 

Bi(x) = ith harmonic amplitude of earth temperature 
cycle at depth x (°C) 

Cj or Cf, a- = Linearized variables (dimensionless) 

D = A least squares estimate of average thermal diffu- 
sivity, a (cm 2 /sec) 

F{6) = Surface temperature function 

d= Amplitude of ith harmonic coefficients = V a 2 -h b 2 t 
(°C) 

M = Total number of temperature data (number of time 
data multiplied by number of depth data) (dimension- 
less) 

Mt = Total number of cycles to be considered (dimen- 
sionless) 

M x = Total number of depths to be considered (dimen- 
sionless) 

N = Highest order of harmonics to be included for 

A expected earth temperature function (dimensionless) 

Pj(x) = Phase angle of ith harmonic of earth tempera- 
ture cycle at depth x, (radians) 

SS= Sum of squares of all of the differences between 
observed and expected earth temperatures 

Si or Sf.A* = Linearized variables (dimensionless) 



T— Period of cyclic data (31,557,600 sec for annual 

cycle) 
t = Earth temperature (°C) 
t x = Earth temperature at x —> °o (°Q 
tk ■ = Observed earth temperature (°C) 
x = Depth in earth (cm) 
xk — Observed depth in earth (cm) 

Greek Characters 

a = Expected thermal diffusivity of earth (cm 2 /sec) 



Bi = Temperature decay factor = \/— i, (cm -1 ) 

al 



/3i = Estimate of fti = \lj^ (cm ~ ] 
. dSS 

^ = Time coordinate used in the text, measured from 

January 1 (sec) 
0fr = Observed time with datum point at January 1 (sec) 
<&i = Phase angle of ith harmonic of earth temperature, 

(radians) <J>; = tan -1 — 
at 

A = Root mean squared deviation (defined in the text). 



1. Objective 

The purpose of this paper is to discuss a harmonic 
analysis of cyclic temperatures at selected depths of 
earth. The existing techniques employing a simple 
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harmonic equation [Kusuda, 1965; Penrod, 1960] or 
Fourier analysis [Carson, 1963] have the disadvantage 
of yielding different thermal diffusivities for a given 
set of data, depending upon the choice of data reduc- 
tion method, choice of cyclic period, choice of two 
depths, and the choice of computational method. The 
technique introduced in this paper eliminates this 
disadvantage and it is statistically more sound than 
previous methods for the analysis of periodic depth 
temperatures of earth. Furthermore, the method 
derived herein may be applied for the accurate calcula- 
tion of thermal diffusivity of the surface layer of other 
celestial bodies in the event cyclic patterns of their 
depth temperatures are obtained in the future. 



2. Introduction and General Background 

Numerous investigators in the past have studied 
earth temperatures with reference to the earth's ther- 
mal diffusivity. A recent publication of Kusuda and 
Achenback [1965] includes an extensive bibliography 
on the subject. Most of the analyses in the papers 
listed in this bibliography apply the well-known peri- 
odic heat conduction formula for either diurnal or 
annual (or both) cycles of observed earth temperatures. 

These analyses have been made to establish the 
average, the amplitude, and the phase angle of cyclic 
temperature patterns at various depth levels, so that 
earth temperatures may be predicted by a computa- 
tional method. The amplitude decays and phase 
angle shifts of temperature cycles with respect to 
depths have been used to estimate the thermal dif- 
fusivity of the station where the earth observations 
were made. The analysis of periodic earth tempera- 
tures, therefore, eventually ends up with four basic 
parameters: thermal diffusivity, average temperature, 
amplitude, and phase angle of the surface temperature 
cycle from a given datum time point. 

The investigators have encountered the following 
problems during the course of their studies. (1) The 
earth surface was difficult to define and earth surface 
temperatures were subject to disturbances from cli- 
matic changes above the surface, such as rainfall, 
cloud cover, evaporation of water, and wind velocity. 
(2) The actual earth temperatures observed did not 
follow an ideal simple-harmonic pattern. Moreover, 
their cyclic patterns varied from year to year and from 
day to day. Except in the work of Carson and in that 
of Kusuda and Achenbach, some degree of arbitrary 
human judgment was exercised to establish meaning- 
ful periodic patterns. For example, when several 
years of monthly average earth temperature records 
were available, the annual thermal diffusivity was 
calculated for each set of annual data. Such calcu- 
lations imply an annual change of thermal diffusivity, 
which is a misapplication of basic heat conduction 
theory, as explained later. (3) For a given station, 
thermal diffusivities were computed both from the 
difference of logarithmic amplitudes and from the 
difference of phase angles for two different depths. 



Unfortunately, the diffusivities calculated by these 
two methods usually did not agree, differing by a 
factor of two in some instances [Kusuda, 1965]. (4) 
When the thermal diffusivities were computed at two 
different pairs of two-depth levels at a given tempera- 
ture station, two different estimates of the diffusivity 
were obtained, even when a consistent method of 
diffusivity computation was employed. 

These four problems resulted in uncertainties of 
thermal diffusivities due to time, due to depth, due to 
method of computation, and due to the method of 
data reduction for a given set of earth conditions at a 
given locality. 

Most of the investigators also attempted to predict 
earth temperatures based on an average of thermal 
diffusivities thus computed, and on other basic param- 
eters; namely, the earth surface temperature ampli- 
tude, the phase angle of the earth surface temperature 
cycle measured from a certain reference datum point, 
and the average earth temperature during the cycle. 

Mathematical description of an earth temperature 
cycle by a Fourier series technique, similar in nature 
with the least squares technique, has been tried by 
several authors. The most recent and comprehensive 
description was made by Carson [1963] for the soil 
of Lemont, 111. The Fourier series technique is basic- 
ally applicable only to a single complete cycle where 
data points consist of a unique set of temperatures 
for equally incremented time coordinates. Carson 
applied the technique to each individual one-year set 
of monthly average earth temperature cycles over a 
period of three years, resulting in three sets of Fourier 
coefficients, which did not agree with each other. 

Kusuda and Achenbach [1965] applied a least 
squares technique in fitting observed earth tempera- 
ture data to simple harmonic expressions for the 
purpose of computing the four basic parameters of 
earth temperature cycles. 

In their least squares technique, Kusuda and Achen- 
bach [1965] incorporated several years' records simul- 
taneously, thus treating the temperature data as 
statistical variables or variates rather than functional 
variables, which is the case in a Fourier analysis. The 
least squares technique fits a sinusoidal curve to the 
observed data in such a manner that the variance 
calculated from the residual deviations of the observed 
temperature from the fitted curve is a minimum. The 
number of data points, the time intervals for the tem- 
perature observations, and the number of data for a 
given time are not restricted as they are in a Fourier 
analysis. 

Kusuda and Achenbach [1965], however, began with 
the premise that the earth temperature cycles follow 
a simple harmonic mode. This premise may be 
reasonably applied to annual cycles of monthly average 
earth temperatures, but it is not quite valid for diurnal 
temperature cycles. They, however, still obtained 
two thermal diffusivities based upon ^he amplitude 
decay and phase angle shift techniques, same as those 
used by other investigators. 

The extension of the least squares technique to cover 
the higher harmonic terms is also desirable if the 
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method is to be applicable to diurnal earth temperature 
cycles of the earth theoretically differing from the 
simple harmonic pattern because of a complex non- 
linear surface heat exchange with the sun and the 
surrounding space. 

In the following section the basic heat conduction 
relation of the periodic temperature cycle is first re- 
viewed. The least squares technique is then applied 
to annual cycles of monthly average earth temperatures 
as an illustrative example. 

3. Basic Heat Conduction Theory on a Semi- 
Infinite Solid 

The surface layer of earth may be assumed to be a 
semi-infinite solid having a homogeneous thermal 
property, as far as the heat conduction theory is con- 
cerned. This assumption, however, may be inade- 
quate in some cases, and a more refined model, such 
as composite solids (an upper layer of finite thickness 
over an infinitely thick solid bed of different thermal 
properties), is desirable. As a theory, this composite 
model probably represents more closely the actual 
surface of earth than the homogeneous model. But 
the thickness of the top layer, and the variation of 
properties between the top layer and the bottom bed 
are usually unknown, which makes the composite 
layer nodel less meaningful than the homogeneous 
model. Also, the mathematical complexity involved 
for the composite solid model may not be justified in 
finding an accurate thermal diffusivity of a surface 
layer from the observed depth temperature cycles. 
Hence, the homogeneous semi-infinite solid model is 
used to calculate average thermal diffusivity repre- 
sented by the surface region of the earth. Or, con- 
versely, temperature cycles at various depths are 
approximated by the equation derived for semi-infinite 
heat conduction models of homogeneous thermal 
properties. 

The heat conduction of the semi-infinite homoge- 
neous solid can be mathematically expressed as 
follows: 



a 



dx 2 



at x 



" dd 
t- 



t a 



atx = t = F{6) 
and F(0+T) = F(6) 



(1) 



(2) 



In actuality it is difficult to define a unique F(0) 
that represents daily or annual cyclic observations of 
earth surface temperatures, since the character of 
these cycles varies from one period to the next. This 
variation is due to the varying nature of heat exchange 
over the earth surface with respect to climatic condi- 
tions. It is postulated that other celestial bodies 
without atmosphere may experience a much more 
consistent and stable F(0) than the earth experiences. 
In this paper it is assumed that there exists a unique 
F(0), which may be a statistical average of many 



periodic cycles. Since F (0) is a periodic function of 
period 7\ it is possible to expand it as follows: 

F(0) = ao+J a, cos (?p OJ + ^bi sin ^f- d\ (3) 

The values of ao, ai, 02, . . . &i, 6 2 , 63 are dependent 
upon the actual nature of the assumed surface tem- 
perature cyclic pattern. 

The cyclic earth temperature [Eckert, 1959] at depth 
x that satisfies (1), (2), and (3) is 



t = t : 



: + J) «>'< 



^cos 



i= 1 



2m0 



sin 



27710 



7^1 



(It should be noted from eqs (3) and (4) that Ao=t x ), 
Formula (4) can also be expressed as 



i= 1 



aT 



(2irid 

—f $l— \\—7^X 



(5) 



where 



Gi=Va 2 + b 2 

tan <$>j = — 

CM. 

By letting 

Gie-yl^j,x=Bi(x)= depth amplitude 



(6) 



6 



(7) 



4>r 



- 1 x — Pi(x) = depth phase angle, 



the thermal diffusivity a can be readily calculated by 
knowing Bj(x) and Pi(x) at two different depths x x and 
x? by the following formulas: 



a hi 



ITT 



, Bi(x 2 ) 

log TV 



BiixA 



£ = 1,2 



(8) 



OLpi = 



17T 

T 



Xx 



Piixi) -Pi(x 2 y 



1, 2 



(9) 



In formulas (8) and (9), notation a B and a P are used 
instead of a to signify that these two thermal diffusi- 
vities are computed by amplitudes and by phase angles 
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of ith harmonic constants, respectively. Theoreti- 
cally, of course 



The best unbiased estimates of a , a u . . . a^, 
61, b 2 . . . 6,v, and Z), that of a are given by the solu- 
tions of the following equations 



atet = api=a for i = l, 2, 3, . . . . 



(10) 



Unfortunately, however, these relations have seldom 
been found in practice. Most previous researchers 
have made extensive studies on am and am only, ig- 
noring higher harmonics. Even for the first harmonic 
constants, agreement between am and an was poor. 
Not only that, the values of a B \ and ap\ varied, as 
mentioned before, depending upon the choice of %\ and 
x 2 of formulas (8) and (9). These difficulties have been 
chiefly attributed to the fact that the earth surface is 
not quite homogeneous and the heat conduction char- 
acteristics and surface temperature cycle of the earth 
surface region vary with respect to physical condition 
and time. The previous difficulties are, however, 
mostly due to the misuse of formula (4). It is not 
correct to use formula (4) either for analysis of one 
given periodic cycle of observation, or for analysis 
of depth variation of thermal diffusitivity. Instead, 
formula (4) should be used for the entire temperature 
data, covering the entire depth and the entire time 
period. For instance, if ten years of monthly average 
earth temperature cycles were complete with five 
depth readings, all of the 600 data points should be 
considered simultaneously for formula (4) to yield 
the best fit values of £«,, (ai, t = 1, 2, 3, . . .) and (bt i i = 1, 2 , 
3 . . .) as well as the best fit value of a. 

It is important here to mention that the foregoing 
statements are based upon the following statistical By letting 
hypothesis: 

Time data (including all the cyclic periods) follow 
basically a given harmonic function of a period T 
(8766 hr for an annual earth cycle). Any deviation 
of actual cycles from a prescribed harmonic function 
are statistical rather than functional. 



dSS 
dAo 

dSS 
dAf 

dSS 



= 



= 0_ 



= (r 



for i= 1,2 . 



N 



(11) 



8- dSS -0 



(11-a) 



where Aq, A\, A 2 , . . . A N , B u B 2 , . . . #v and D are 
the best estimates of ao, ai, a 2 , • . . a#, 60, 61, b 2 , . . . 
6a, and a. 



M T N /2tt7 - \ 

And SS=2 k-Ao-^A^kcoB^ek-fitt) 



-g^e-^sin^^,-^ 



N 

2 

1=1 



(12) 



vhere 



1-4 



ITT 

DT 



Ci, k - = e~^i x k cos f— Ok — PiXk) 



(13) 
(14) 



4. The Least Squares Technique 

Although the technique described in this paper 
applies to any kind of periodic temperature data ob- 
served at several depths of any semi-solid system, 
the specific treatment will be given for annual cycles 
comprising monthly average earth temperatures taken 
at selected depths. 

Assume that the number of annual cycles of earth 
temperature is M* and each annual cycle consists of 
twelve monthly average readings at each of M x depths. 
If there are no missing data, although this is not a 
requirement, the total number of observations M 
should be 

M=M t xM x xl2. 

Equation (4), however, should be modified to include 
only the finite number of harmonic terms up to some 
larger number N. In other words, it is assumed that 
Oi and hi will be numerically insignificant for all i > N. 



and by noting that 
dSS 



dSS 



is also satisfied by rr— =0 as long as ZM 0, 
ou dpi 

the equations in (11) can be expressed simply by the 
linear matrix equation shown in table 1. The left- 
hand side of the equal sign of table 1 is actually the 
product of two matrices, one square matrix of (27V 
-f 1) X (2/V-f 1) and one column matrix of dimension 
(2/V-f-l). Equation (ll~a) is nonlinear and is ex- 
pressed in terms of linearized variables, C^a and S;,a 
as is shown in table 2. 

Determination of Ao, A\, A 2 , . . . An, B\, B 2 , . . . 
B N as a best estimate of a , ai, a 2 , . . . a^, 61, b 2 , . . . 
by is straightforward if the (w,a's and S/,a-'s are known. 
However, these variables include an unknown D which 
is to be solved by (ll~a). 

The method employed in this paper is to calculate 
a number of sets of /4 , Au B\\ i= 1, 2, . . . /V, for sev- 
eral assumed D's. These calculated sets are then 
applied to the equation in table 2. By an iterative 
procedure, a set of A , Ai, J5,; i= 1, 2, . . . N, and D 
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Table 1. Matrix expression of least squares equations (normal equations) 



Least Squares Equations 



M M 

j~l A = l 

(ixl) (1 x/V) (1 x/V) 


M M M 

]T Ci t k 2 Ci,kCj y k ^ Ci, k Sj,k 

A = l A = l A-l 

(N xi) (NX N) (TV x N) 

M M M 

2 Si.fr ^ Si t k-Cjk 2 &, A'Sj, A- 

A = l A = l A = l 

(/Vxl) (/Vx/V) (/VXiV) 
i 



X 



,4 



(ixl) 




B, 



(/Vx 1) 



Ci,/,=e'yj~n cos (-=- ft. — "W 7™ 



Si, a = e- V^. >« sin f -=- ft — VTJ^ **) 



I 

A = l 



2> 



(1X1) 



2 c«.*«* 



(A' xi) 



^ Si, ktk 



(NX I) 



TABLE 2. Mathematical formula for the partial derivative of the sum of squares of all of the differences between observed and expected earth 

temperatures with respect to thermal diffusivity 



2 VKAi + Bi) J «***Ct.*-2; VJiA-Bi) J fc**S,, ft 

' A- / A- 

' A i A- 

JV V M A V ^ M 

2 ^ 2 vj(4+Bj) 2 Ci, k c J , k x k + 2 4 2 V/U-fl;) 2 c *. *«>.*** 

' J A' i j A- 

JV JV A/ tf W ^_ M 

2 fii 2 V/C4+*/) 2 s *. *Q. *** +2 Bi I vj^-iy 2 Si, *s Jf *** 



= 6-0 

can be found to satisfy 8 = (see table 2). These co- The value of A, or the square root of the average of 

efficients and the assumed D are then substituted the sum of squares of deviations is called hereafter 

into formula (12) for the calculation of the mean the root mean squared deviations (RMSD) of the ob- 

squared deviation of the fitted least square values from served temperatures from those calculated by eq (4). 

the observed data by the following equation: The term "standard deviation"' is purposely avoided 

here since the properties of these deviations resulting 

^ 2 __oo from the iterated nonlinear fitting procedure are not 



M 



known. 
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Table 3. Observed earth temperature data of Lemont, 
Carson 



III., by 











Depth (cm) 








Month 
































1 


10 


20 


50 


100 


305 


884 


1 


-0.4 


-0.2 


0.2 


1.8 


4.0 


9.7 


11.3 


2 


-0.2 


-0.3 


0. 


1.2 


3.0 


8.3 


1 1 .3 


3 


4.0 


3.3 


3.0 


2.9 


3.5 


7.3 


11.4 


4 


8.1 


7.1 


6.8 


6.3 


6.0 


7.2 




5 


16.2 


15.5 


14.4 


12.4 


10.2 


7.9 




6 


22.4 


22.5 


21.3 


18.7 


15.5 


9.7 




7 


24.6 


24.4 


23.5 


21.4 


18.7 


12.0 


10.1 


8 


24.0 


23.7 


23.0 


21.7 


19.8 


13.8 


10.1 


9 


19.2 


18.8 


19.0 


19.1 


18.8 


14.9 


10.3 


10 


14.1 


14.2 


14.6 


15.3 


15.9 


14.8 


10.4 


11 


6.9 


7.2 


8.2 


10.1 


12.1 


13.9 


10.6 


12 


1.8 


1.9 


2.8 


4.8 


7.8 


12.5 


11.0 


13 


-0.6 


-0.5 


0.2 


1.8 


4.6 


10.5 


11.1 


14 


1.1 


0.8 


0.9 


1.8 


3.5 


8.9 


11.2 


15 


2.4 


2.1 


2.3 


2.9 


3.9 


8.0 


11.2 


16 


10.7 


10.1 


9.3 


7.9 


6.7 


7.5 


11.0 


17 


14.9 


14.2 


13.4 


12.1 


10.6 


8.4 


10.7 


18 


23.6 


22.9 


21.4 


18.6 


15.4 


10.0 


10.6 


19 


25.3 


24.7 


23.8 


21.7 


19.0 


12.3 


10.4 


20 


23.9 


23.2 


22.8 


21.4 


19.7 


14.2 


10.5 


21 


20.0 


19.7 


19.7 


19.4 


18.9 


15.1 


10.6 


22 


13.5 


13.3 


14.0 


15.2 


16.1 


15.2 


10.9 


23 


6.2 


6.0 


6.9 


8.7 


10.9 


14.2 


11.1 


24 


1.3 


1.1 


1.9 


4.0 


7.2 


12.3 


11.3 


25 


0.5 


0.4 


1.0 


2.7 


4.8 


10.3 


11.4 


26 


0. 


-0.4 


-0.1 


1.2 


3.4 


9.0 


11.5 


27 


2.9 


2.1 


2.1 


2.6 


3.6 


7.8 


11.5 


28 


11.3 


11.0 


10.2 


8.6 


7.1 


7.4 


11.1 


29 


16.7 


16.5 


15.7 


13.9 


11.8 


8.7 


11.0 


30 


20.6 


20.1 


19.3 


17.4 


15.0 


10.4 


10.7 


31 


26.7 


26.0 


24.7 


21.8 


18.7 


12.3 


10.5 


32 


25.5 


25.1 


24.4 


22.9 


20.7 


14.2 


10.5 


33 


19.5 


19.2 


19.3 


19.3 


19.0 


15.3 


10.7 


34 


13.3 


13.1 


13.6 


14.7 


15.7 


15.3 


10.8 


35 


4.5 


4.3 


5.2 


7.7 


10.6 


14.0 


11.0 


36 


-0.3 


-0.4 


0.4 


2.7 


5.8 


12.1 


11.4 



Sample Calculations: The normal equations of 
tables 1 and 2 have been solved for annual earth tem- 
perature data of Lexington, Ky., and Lemont, 111., 
which had been obtained by Penrod [1960] and Carson 
[1963], respectively. The Kentucky data covered five 
years, from 1952 through 1956, at six depths, 0, 61, 
122, 183, 244, and 305 cm. Except for the data at 61 
cm and 122 cm during July and August of 1953, the 
data are complete. Lemont data shown in table 3 are 
complete for the period of three years, from 1953 
through 1955, at seven depths, 1, 10, 20, 50, 100, 305, 
and 884 cm, except at the depth of 884 cm during May, 
June, and July of 1953. 

Figures 1 and 2 illustrate the trend of 8 (table 2) 
and corresponding values of A calculated for several 
values of thermal diffusivities for the earth temperature 
data of Carson and Penrod, respectively. The calcu- 
lations in figures 1 and 2 were carried out for TV =4. 

It is seen in figure 2 that the value of A, the RMSD 
of calculated temperatures from the observed data, 
decreases rapidly to the value of 1.07 °C as the assumed 
thermal diffusivity increases to a value at 0.005 
cm 2 /sec. On the other hand it is seen from figure 2 
that the diffusivity values that satisfy 8 = for the 
Lemont data is a = 0.0058 cm 2 /sec. Also it will be 
interesting to compare these thermal diffusivities with 
those calculated by conventional techniques by means 
of amplitude decay and phase angle shift by Penrod 
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FIGURE 1. ^ = ~^o~ and A plotted against thermal diffusivity for 
earth temperature data of Lemont, III. 
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FIGURE 2. & = -fijf and A plotted against thermal diffusivity for 
earth temperature data of Lexington, Ky. 



48 



from his Lexington data and with the data deduced 
from Carson's data for Lemont, III. 

Lexington data a ln = 0.0088 cm 2 /sec 
an =0.0056 cm 2 /sec 

Lemont data a/ n = 0.0048cm 2 /sec 
a P1 =0.0045cm 2 /sec. 

Although the thermal diffusivities calculated by con- 
ventional amplitude and phase angle techniques differ 
from one another considerably, the values of A cor- 
responding to their diffusivities are all less than 1.3 °C 
in their ranges, as seen in figures 1 and 2. (This small 
change of A with respect to a wide range of thermal 
diffusivity may be the very reason that the two methods 
can yield quite different thermal diffusivities.) 

The mean squared deviation A 2 may not be as in- 
sensitive to the thermal diffusivity for the other types 
of periodic temperature data as it has been found to 
be in the analysis on annual cycles of monthly average 
earth temperatures. 

The Lexington earth temperature data were used 
to examine the influence of including higher harmonic 
terms on the values of A and it was found that the 
least squares analysis, including the higher harmonic 
terms up to degrees 2, 4, and 8, do not yield values 
of A which are significantly different from the results 
of the analysis ignoring all of the higher harmonics 
for all the thermal diffusivities assumed. 

Sample sets of least squares constants calculated 
for these two stations are shown in table 4. 



TABLE 4. Least squares constants for earth temperature data of 
Lexington, Ky., and of Lemont, III. 





Lexington. K\. 


1 ,emont, III. 


M 


346 


249 


N 


4 


4 


1) 


0.0050 


0.0058 


A 


1 .05 


.72 


/„ 


14.48 


1 1 .39 


•/, 


- 10.58 


- 10.63 


A, 


-O.30 


0.08 


A 3 


- .54 


-.16 


A 4 


.24 


-.30 


Hi 


-7.13 


-7.97 


Bt 


-0.38 


0.59 


B 3 


0.37 


.21 


B 4 


.14 


.39 



It may be interesting to compare the calculated earth 
temperatures based upon the least squares constants 
with the original observed data. Table 5 shows the 
calculated earth temperatures when N=4> is employed 
for formulas of table 1 on the observed data of Lemont. 
111. (table 3). Similar calculations such as those for 
table 4 for /V=l were also made to evaluate the im- 
provement of taking higher harmonics over using the 
simple harmonic expression for the annual cycle of 
Lemont data. According to such comparison changes 
in the calculated results due to the inclusion of higher 



harmonic terms are very small. Table 5 also may be 
compared with table 6, which is the arithmetic average 
(or norm) of the 3-year cyclic data of table 3. Sig- 
nificant discrepancy of the arithmetic average from the 
least squares result exists for the region near the 
surface. 



TABLE 5. Predicted earth temperatures for Lemont. Illinois, using 
least squares constants up to fourth harmonics 











Depth (cm) 








Month 
































1 


10 


20 


50 


100 


305 


881 


1 


-0.5 


-0.1 


0.4 


1.9 


4.1 


10.5 


11.7 


2 


-0.3 


-0.0 


0.3 


1.4 


3.1 


9.0 


11.7 


3 


2.6 


2.6 


2.6 


2.8 


3.6 


8.0 


11.7 


4 


9.4 


9.0 


8.6 


7.6 


6.7 


7.7 


11.6 


5 


15.9 


15.3 


14.7 


13.0 


10.9 


8.5 


11.4 


6 


21.8 


21.1 


20.3 


18.2 


15.4 


10.0 


11.3 


7 


25.3 


24.7 


24.0 


22.1 


19.2 


12.0 


11.1 


8 


23.7 


23.5 


23.2 


22.3 


20.5 


13.9 


11.0 


9 


19.2 


19.2 


19.3 


19.4 


19.0 


15.0 


11.1 


10 


13.2 


13.6 


14.0 


15.0 


15.9 


15.1 


11.2 


II 


5.3 


6.0 


6.7 


8.7 


11.1 


14.2 


11.4 


12 


0.5 


1.1 


1.8 


3.7 


6.6 


12.5 


11.5 



V/ = 249. 

V = 4. 
#=.0058. 



Table 6. Arithmetic average earth temperature of Lemont , Illinois. 
earth temperature data in table 3 











Depth (cm) 








Month 
































1 


10 


20 


50 


100 


305 


884 


1 


-0.2 


-0.1 


0.5 


2.1 


4.4 


6.8 


11.3 


2 


.3 


.0 


.3 


1.4 


3.3 


8.7 


11.3 


3 


.0 


3.1 


2.5 


2.8 


3.7 


7.1 


11.4 


1 


10.0 


9.4 


8.8 


8.7 


6.6 


7.4 


li.l 


5 


15.0 


15.4 


14.5 


12.8 


10.0 


8.3 


10.8 


6 


22.2 


21.8 


20.7 


18.3 


15.3 


10.0 


10.7 


7 


25.5 


25.1 


24.0 


21.6 


18.8 


12.2 


10.3 


8 


24.4 


21.0 


23. 1 


22.0 


20.1 


i 1.1 


10.4 


9 


19.6 


I'). 2 


10.3 


I'). 3 


18.9 


16.1 


10.6 


10 


13.6 


13.6 


14.1 


15.1 


15.9 


15.1 


10.7 


II 


5.9 


5.8 


6.8 


8.8 


11.2 


14.1 


10.9 


12 


0.9 


0.9 


1.7 


3.8 


6.9 


12.3 


11.2 



Figure 3 compares the least squares calculation of 
the first harmonic constants with those derived from 
the observed three annual cycles of earth temperature 
in Lemont, 111. 

The coordinates a x and b x of the figure represent 
first harmonic constants describing the observed earth 
temperature at a depth x from the ground surface. 
Since the Lemont data consists of seven depth observa- 
tions for a period of three years, each depth point is 
indicated with three numerals representing the first, 
second, and third year. 

The solid curve was obtained by connecting a x and 
b x calculated by theoretical relation (depicted in the 
upper right corner of fig. 3) for several depths. The 
calculations of a x and b x were made with the least 
squares constants indicated in table 5. The curve 
representing the computed values of a x and b x lit 
extremely well for all of seven groups of three points. 
A similar comparison was made for the second har- 
monic constants of the same earth temperature data. 



49 



Calculated by 



a x = e -/3x |- A| cos(-^-,3 x ) + B l sin (i--^)]- 
b x =e"^ x [-A,sin(-J-i8x)+B l cos (-|-^ x )]- 



T=365,25 X24 X 3600 sec 




FIGURE 3. Graphical representation of harmonic constants for the earth temperature cycles 

of Lemont, III. 



Although the scatter of these second harmonic con- 
stants derived from observed data is much more than 
those indicated in figure 3, the computed curve did 
closely follow the centroid of the three points. 

5. Conclusions 

With the use of annual cyclic earth temperature 
data of Lexington, Ky., and of Lemont, 111., a new 
technique is illustrated to yield single thermal diffu- 
sivity that will satisfy the least squares relation from 
the entirety of observed data. The technique requires 
an iterative solution of a nonlinear part of normal 
equations. The calculated earth temperatures based 
upon the theoretical heat conduction relation agrees 
very well with the average annual cyclic data of the 
monthly average earth temperatures when the least 
squares constants are properly evaluated. 

It has also been found in this analysis that the devia- 
tion of the observed earth temperatures from the 
calculated temperatures is fairly insensitive to the 
thermal diffusivity as long as the thermal diffusivity 
is in the neighborhood of 0.005 cm 2 /sec, at least for 
the annual cycles of monthly average earth tempera- 
tures for two stations studied in this paper. This 
insensitivity may be the very reason that the current 
and existing techniques resulted in two quite different 
thermal diffusivities depending upon the method of 
calculation used, the amplitude method or the phase 
angle method. 



The method developed in this paper is recommended 
as the preferred method for the determination of 
thermal diffusivity from the periodic data of depth 
temperatures, not only of the earth surface by of other 
celestial bodies as well, in the event such data becomes 
available. 



This paper is a modification of a report prepared 
under the contracted research for the Office of Civil 
Defense. The permission from OCD to release the 
paper for publication is appreciated. 

The author greatly appreciates much critical review 
given by A. McNish of the National Bureau of 
Standards. 
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